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A two-step unified framework for the evaluation of continuum field expressions from 
molecular simulations for arbitrary interatomic potentials is presented. First, point- 
wise continuum fields are obtained using a generalization of the Irving-Kirkwood 
procedure to arbitrary multi-body potentials. Two ambiguities associated with the 
original Irving-Kirkwood procedure (which was limited to pair potential interactions) 
are addressed in its generalization. The first ambiguity is due to the non-uniqueness 
of the decomposition of the force on an atom as a sum of central forces, which is a re¬ 
sult of the non-uniqueness of the potential energy representation in terms of distances 
between the particles. This is in turn related to the shape space of the system. The 
second ambiguity is due to the non-uniqueness of the energy decomposition between 
particles. The latter can be completely avoided through an alternate derivation for 
the energy balance. It is found that the expressions for the specific internal energy 
and the heat flux obtained through the alternate derivation are quite different from 
the original Irving-Kirkwood procedure and appear to be more physically reasonable. 
Next, in the second step of the unified framework, spatial averaging is applied to the 
pointwise field to obtain the corresponding macroscopic quantities. These lead to 
expressions suitable for computation in molecular dynamics simulations. It is shown 
that the important commonly-used microscopic definitions for the stress tensor and 
heat flux vector are recovered in this process as special cases (generalized to arbitrary 
multi-body potentials). Several numerical experiments are conducted to compare the 
new expression for the specific internal energy with the original one. 
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I. INTRODUCTION 


The idea of defining continuum fields from particle mechanics (for the special case of 
pair potential interactions) was pioneered in the landmark paper of Irving and Kirkwood--. 
Irving and Kirkwood derived the equations of hydrodynamics from the principles of non¬ 
equilibrium classical statistical mechanics and in the process established pointwise definitions 
for various continuum fields. Under this procedure, basic continuum fields including the 
mass density, momentum density and the specific internal energy are defined a priori using 
a probability density function. Using these definitions, expressions for the stress tensor and 
the heat flux vector fields are obtained that identically satisfy the balance laws of continuum 
mechanics. The continuum fields obtained in Irving and Kirkwood’s original paper- involved 
a series expansion of the Dirac delta distribution, which is not mathematically rigorous.- 
In a follow-up study, Noll^ proved two lemmas, which allowed him to avoid the use of the 
delta distribution and to obtain closed-form analytical expressions for the continuum fields. 

Since the Irving-Kirkwood procedure is stochastic in nature, many problems arise when 
one tries to use the resulting expressions for a practical calculation — a key one being our 
lack of knowledge of the probability density function. To avoid these difficulties, Hardy^ 
and independently Murdoch^""- developed a simpler spatial averaging procedure that avoids 
the mathematical complexity of the Irving-Kirkwood procedure. We refer to the procedure 
due to Hardy as the Hardy procedure and that due to Murdoch as the Murdoch procedure .— 
In these procedures, continuum fields are defined as direct spatial averages of the discrete 
equations of motion using a normalized weighting function. This approach also leads to a 
set of definitions that identically satisfy the balance equations. Therefore, we have three 
different approaches for defining the continuum fields from particle mechanics — although 
originally developed for pair potentials only. 

Of the continuum fields, the stress tensor has been studied most extensively. In addition 
to the definitions for the stress tensor obtained from the systematic approaches described 
above, a number of other definitions have been proposed in the past dating back to the work 
of Cauchy— on the stress vector and Clausius^ on the virial stress.— Efforts at obtaining 
microscopic definitions for the stress tensor (as well as other continuum variables) are on¬ 
going; see for example Refs^ - — for some important contributions. A recent article^ by the 
authors extensively studies the definition for the stress tensor within a unified framework 
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based on a generalization of the Irving-Kirkwood procedure to arbitrary multi-body poten¬ 
tials followed by a process of spatial averaging. Through this unified framework it is shown 
that all existing definitions, including the virial stress tensor,— Hardy stress tensor,- and 
Cauchy/Tsai stress tensor, 11,12 ' 15 ' which all seem to be derived from disparate approaches, 
follow as special cases from a single stress expression. Furthermore, the derivation in Ref 2 ' 
reveals the subtle (and hitherto unrecognized issue) that interatomic potentials constitute 
continuously differentiable extensions to functions defined over a more limited domain. This 
is a vital part of the derivation with important implications for the uniqueness of the mi¬ 
croscopic stress tensor — an issue which is widely discussed in the literature cited above. 
Although there have been a number of attempts to generalize the Irving-Kirkwood proce¬ 
dure and the Hardy procedure to multi-body potentials (see Refs^ 2 ’— J2), these attempts 
are either restricted to specific potentials (see Refs^^) or the source of non-uniqueness 
of the stress tensor is not explicitly identified. In contrast, the unified framework devel¬ 
oped in Ref 27 applies to arbitrary multi-body potentials and rigorously characterizes the 
non-uniqueness of the stress tensor. 

The aim of this paper is to continue to use this unified framework to study the energy 
balance equation of continuum mechanics in the context of multi-body potentials. As noted 
earlier, in the original Irving-Kirkwood and the Hardy procedure, the definition for the 
potential part of the specific internal energy (for the special case of pair potentials in a 
mono-atomic system) is assumed a priori and the expression for the heat flux vector is then 
derived to ensure that the energy balance equation is identically satisfied. Unfortunately, 
this approach does not generalize to arbitrary multi-body potentials (or even pair potentials 
with multiple species types) since it involves an ambiguous definition for the “energy of 
an atom”. To the best of the authors’ knowledge, all the existing works (see Refs^^—) 
which attempt to derive a microscopic definition for the heat flux in the case of multi-body 
potentials by generalizing the Irving-Kirkwood procedure or the Hardy procedure suffer 
from this ambiguity. For example, in Ref^ 2 it was assumed that the energy corresponding to 
a cluster of three particles interacting through a three-body potential is evenly distributed 
among the particles. However, there is no symmetry argument to justify this assumption.— 
Furthermore, even for the case of pair potential interactions, the original Irving-Kirkwood 
approach leads to an expression for the heat flux vector which is not invariant with respect 
to the addition of a constant to the potential energy of the system, which is not physically 
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reasonable. In contrast, in the Murdoch procedure, the specific internal energy and heat 
flux vector are obtained together as part of the derivation and the resulting expressions are 
consistent with physical expectations. Motivated by this, in this paper, we reformulate the 
Irving-Kirkwood procedure using the method followed by Murdoch^. This approach leads to 
physically-acceptable expressions for the internal energy density and heat flux vector which 
are grounded in rigorous statistical mechanics principles and which does not require any 
energy decomposition between the particles. Furthermore, as noted above, our derivation 
extends those of Irving-Kirkwood and Murdoch to arbitrary multi-body potentials. Finally, 
the application of the spatial averaging step in the unified procedure leads to expressions 
suitable for use in molecular dynamics simulations. These expressions are compared with 
those from the original Irving-Kirkwood formulation through a number of simple numerical 
experiments. 

The following notation is used in this paper. Vectors are denoted by lower case letters in 
bold font, while tensors of higher order are denoted by capital letters in bold font. The inner 
product of two vectors is given by a dot , and their tensor product is given by the symbol 
“(g)”. The inner product of two second-order tensors is denoted by The gradient of a 
vector field, v(x), is denoted by V x v(x). A second-order tensor, T, operating on a vector, 
v, is denoted by Tv. The divergence of a tensor field, T(x), is denoted by div* T(x), which 
corresponds to dTij/dxj in indicia! notation (with Einstein’s summation convention). 


II. CONTINUUM FIELDS AS PHASE AVERAGES 

Consider a system modeled as a collection of N point particles, each particle identified by 
an index a (a = 1, 2,..., N). The position, mass, and velocity of particle a are denoted by 
x a , m a and v Q , respectively. We assume that the particles interact through a continuously 
differentiable function V{x \,..., £C/v), which is called the potential energy of the system. The 
complete microscopic state of the system at any instant of time is known from the knowledge 
of position and velocity of each particle in M 3 . Hence, the state of the system at time t may 
be represented by a point in a 6iV-dimensional phase space.— Let T denote the phase space. 
Therefore any point in F, can be represented as, 

(x(t); v(t)) := (xi(t ),..., x N (t)] v x (t ),..., v N (t)). (1) 
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In reality, the microscopic state of the system is never known to us and the only observables 
identified are the macroscopic fields as defined in continuum mechanics. We identify the con¬ 
tinuum fields with macroscopic observables obtained in a two-step process: (1) a pointwise 
held is obtained as a statistical mechanics phase average; (2) a macroscopic held is obtained 
as a spatial average over the pointwise held. The phase averaging in step (1) is done with 
respect to a continuously differentiable^ probability density function W : T x M + —>- M + 
defined on all phase space for all t. The explicit dependence of W on time t, indicates that 
our system need not be in thermodynamic equilibrium. 

The basic idea behind the original Irving and Kirkwood procedure is to prescribe the 
mass density, velocity and the specihc internal energy helds, which we call the input helds, 
and derive the body force vector, stress tensor and the heat hux vector helds, which we call 
the output helds, such that all the definitions are consistent with the balance laws of mass, 
momentum and energy: 


Input fields 


Output fields 


\ 

mass density 


/ 

body force 

velocity 

> ->■ < 

stress 

specific internal energy 


heat hux 


( 2 ) 


To arrive at these definitions, we repeatedly use the following result of Liouvillc’s theorem, 
which describes the evolution of the probability density function: 

N 


dW 


+ 'y ^ \ v a • v a , ct iK + v a • v„ a ir ] — o. 


( 3 ) 


a=l 


Since the force on a particle a is given by 

fa ■= -V^V, 


( 4 ) 


equation (J3]) can be rewritten as 

N r 


dW 

dt 


E 

a=l L 


v, v 

v a -V x W-^^-V v W 


nir 


= 0, 


( 5 ) 


where, as stated before, V(a? 1; x 2 ,..., x N ) denotes the potential energy of the system. Equa¬ 
tion ((5]) is called Liouville’s equation. 

To proceed, we divide the potential energy into two parts: 
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1. An external part, V ex t, associated with long-range interactions such as gravity or 
electromagnetic fields, 

2. An internal part, Vi n t, associated with short-range particle interactions. In general, 
the internal part of the potential energy is also called the interatomic potential energy. 

We next define the input fields used in the Irving-Kirkwood procedure. 


A. Phase averaging 


Under the Irving-Kirkwood procedure, pointwise fields are defined as phase averages. 
For example, the pointwise mass density held is defined as 


p(x,t) := Vm a 

Jr 3 


W8{x a — x ) dxdv, 


( 6 ) 


Ir 3 N xR 3N 

8 denotes the Dirac delta distribution, and denotes summation from a = 1 to N. To 
avoid the Dirac delta distribution and for greater clarity we adopt the notation introduced 
by Noll. Hence (JHJ) can be rewritten as 

p{x,t) = E m a / W dx i... dx a _idx a+ i... dx N dv 

a ^ 

=: ^m Q (W \ x a = x), (7) 

a 

where ( W \ x a = x) denotes an integral of W over all its arguments except x a , and x a is 
substituted with x. 

The second input held, which is the pointwise velocity held, is dehned via the momentum 
density held, p(x,t), as follows: 

p(x,t) \= (Wv Q | x a = x) , (8) 

a 

p(x,t) 


v(x,t) : = 


( 9 ) 


p(x,t) 

The third input held, which is the specihc internal energy, depends on the interatomic 
potential. At this point, it must be noted that the original Irving-Kirkwood procedure was 
limited to systems interacting through a pair potential function: 


Vmt — Vi n t(ri2, - - -, tin , r 23 ,..., r-(jv—tjjv) 

= E v «> 


( 10 ) 
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where V a is the energy of particle a, defined as 


V Q := 


1 

2 


^ ^ 4 > /3a{ r l3a) + ^ ^ 4 > ap{ r ap) 

- P P 

/3<a P>a 


( 11 ) 


and 0 Q( g (a < /3) is the energy corresponding to the interaction of the pair (a,/3). In this 
case, the specific internal energy is defined as 

e(x,t) ■= e k (x,t) + e v (x,t), (12) 


where 


pck(x,t) 


^2m a (\\v a \\ 2 W | x a = x), 


is the kinetic contribution to the specific internal energy , and 


(13) 


pe v (x, t) - x,, = x), (14) 

a 

is the potential contribution to the specific internal energy. According to the definition given 
in (fT2]h the specific internal energy at (x, t) is the weighted sum of the energy of each particle 
with the probability that it is at x at time t. It is clear from the definition in (Hill that the 
interaction energy (j) a p, between any two particles a and /3, is shared equally between the 
particles a and (3. This is plausible for systems with identical particles interacting with pair 
potential, but there is no a priori physically motivated way of deciding how to distribute the 
energy for systems interacting through a multi-body potential. This is one of the primary 
reasons why the definition for the specific internal energy and the energy balance equation 
has to be re-examined as we do later in Section III D1 

It is clear from the definitions in (J7J) , (JEJ), (fT3l) and (fT4l) that the integrals in these 
equations converge only under appropriate decay conditions^! on W. Under these condition, 
any continuously differentiable vector or tensor-valued function defined on the phase space 
for all t (and satisfying certain additional decay conditions described in Ref 27 ), we have^ 


^Xa fU dx a 

/ W div a , a Gdx a , 

(15a) 


J R3 


V Va W dv a = - 

/ W div„ a G dv a . 

(15b) 


The above identities are repeatedly used in deriving the equation of continuity, the equation 
of motion, and the energy balance equation in the Irving-Kirkwood procedure. 
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B. General interatomic potentials 


In this section, we describe some properties of interatomic potentials, which play a cru¬ 
cial role in extending the original Irving-Kirkwood procedure to multi-body potentials. In 
addition, it gives new new insights into the original procedure which was limited to pair 
potentials. This section is largely based on Ref 2 -, and is briefly summarized here for com¬ 
pleteness and to define the necessary notation and terminology. 

In general, the internal part of the potential energy, also called the interatomic potential 
energy, depends on the positions of all particles in the system: 

k^int ldnt(^T; *®2, • • • , ^Tv), (16) 

where the “hat” indicates that the functional dependence is on absolute particle positions 
(as opposed to distances later on). We assume that Vj nt : M 3JV —> R is a continuously 
differentiable function.— Due to the invariance of the potential energy with respect to rigid- 
body motions and reflections, it can be shown that Vi nt in can be expressed as a new 
function 37 

V int = Vint(-), (17) 

where the argument of V int is an N(N — l)/2 tuple of “physically-realizable distances”. 
Before we describe what this means, we note that the N(N — l)/2 distances between the 
N particles embedded in M 3 are not independent. This can be easily seen for any collection 
of 5 particles or more. Therefore, the set of all N(N — l)/2 tuples of physically realizable 
distances is a proper subset of ]jW(iV-i)/ 2 . In fact, it is a (3 N — 6)-dimensional manifold 
called the shape space of the system which is defined as 

S '■= {( r i2, ri 3 ,... ,riN,r 2 3, ■ ■ ■, ^(jv-i)at) I 

r a p = \\x a -xg\\, (®i,..., x N ) e R 3N }. (18) 

For example, for a chain of 3 particles in one dimension, with positions X\ < x 2 < x 3 , the 
three distances (^ 12 ,^ 13 , ^ 23 ) must satisfy r 12 +r 23 = r 13 . Values of (^ 12 ,^ 13 ,^ 23 ) that do not 
satisfy this constraint are not “physically realizable” and are therefore outside the shape 
space manifold. 

From the above discussion it is clear that the potential energy is only defined on the 
shape space of the system. We will soon see that in order to derive the stress tensor, we 
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need to evaluate partial derivatives like the following: 


dr 12 


lim 

e —>0 


Mnt (r 12 + e,... ,r jv(jv— i) ) — Vi n t(ri 2 , ■ ■ ■, r n(n- i)) 


( 19 ) 


e 


It is clear that this relation requires us to evaluate the potential energy outside the shape 
space since if (r i2 ,..., rjvpv-i)) is on S then by adding e to one of the distances, we move 
off it. Thus, the expression in (TT9|) makes sense only when we extend the function to the 
neighborhood of the shape space manifold (see Section 3.4 in Ref2£ for a more detailed 
discussion). 

This is the reason we now restrict our discussion to those systems for which there exists 
a continuously differentiable extension of V int , defined on the shape space, to M. N ( N ~d/ 2 . 
This is a reasonable assumption because all interatomic potentials used in practice, for 
a system of N particles, are either continuously differentiable functions on ]gA(v-i)/ 2 ^ or 
can easily be extended to one. For example, the pair potential and the embedded-atom 
method (EAM) potential™ are continuously differentiable functions on fl£-W(JV-i)/ 2 ^ while 
the Stillinger-Weber— and the Tersoff— potentials which depend on the angles between 
relative position vectors, can be easily extended to by expressing these angles 

as a function of distances between particles. Therefore, we assume that there exists a 
continuously differentiable function Vi nt : M^^ -1 '/ 2 ]g^ such that the restriction of V; nt to 

S is equal to V int : 


Vi n t(s) = Vmt(a) Vs = (r 12 ,..., r (iV - i)n) e S. 


( 20 ) 


An immediate question that arises is whether this extension is unique in a neighborhood 
of s G S. Note that for 2 < iV < 4, 3N — 6 = N(N — 1)/2. Therefore, for 2 < iV < 4, for 


every point s G S, there exists a neighborhood in W. N ^ N d/ 2 which lies in S. However, for 
iV > 4, there may be multiple extensions of V int . 


We will soon see that the quantity evaluated in (|T9|) may differ for different extensions. 
On the other hand, the internal force on any particle a, 



( 21 ) 
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is uniquely defined for any extension. We next address the possibility of having multiple 
extensions for the potential energy by studying the various constraints that the distances 
between particles must satisfy in order to be embeddable in M 3 . We demonstrate, through 
a simple example, how multiple extensions for the potential energy lead to a non-unique 
decomposition of the force on a particle, which in turn leads to a non-unique pointwise 
stress tensor. 

Central-force decomposition and the possibility of alternate extensions 

We first show that the force on a particle can always be decomposed as a sum of central 
forces regardless of the nature of the interatomic potential. The force on a particle due to 
internal interactions is defined in (12 1 11 . This can also be evaluated using the continuously 
differentiable extension Vi nt and the chain rule as 


fa\ r 12, • • •, J'(jv-i)jv) = -V 3 , Q V int (r 12 , • • •,r (J v-i)iv) 



( 22 ) 


P 

P ¥«* 


where 



(23) 


is the contribution to the force on particle a due to the presence of particle (3. 

Note that f a p is parallel to the direction xp — x a and satisfies f a p = —fp a - This leads us 
to the important result that the internal force on a particle, for any interatomic potential 
that has a continuously differentiable extension, can always be decomposed as a sum of central 
forces, i.e., forces parallel to directions connecting the particle to its neighbors. This may 
seem strange to some readers due to the common confusion in the literature of using the 
term “central-force model” to refer exclusively to simple pair potentials. In fact, we see that 
due to the invariance requirement stated above, all interatomic potentials (including those 
with explicit bond angle dependence) that can be expressed as a continuously differentiable 
function of distance coordinates, are central-force models. By this we mean that the force 
on any particle (say a) can be decomposed as a sum (over (3) of terms, f a p, aligned with 
the vectors joining particle a with its neighbors and satisfying action and reaction. The 
difference between a pair potential and a many-body potential is that in the former f a y 
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only depends on r a p whereas in the latter f a p can depend on the distances between all 
particles. 

The next question is how different potential energy extensions affect the force decompo¬ 
sition in f|22|) . We have already seen through (12T| that the force /“* is independent of the 
particular extension used. However, we show below that the individual terms in the decom¬ 
position, fap, are not unique. These terms depend on the manner in which the potential 
energy, defined on the shape space, is extended to its neighborhood in K^^ 1 )/ 2 , 

In order to construct different extensions, we use the geometric constraints that the 
distances have to satisfy in order for them to be embeddable in R 3 .— The nature of these 
constraints is studied in the field of distance geometry, which describes the geometry of sets 
of points in terms of the distances between them. One of the main results of this theory, 
is that the constraints are given by Cayley-Menger determinants, which are related to the 
volume of a simplex formed by N points in an N — 1 dimensional space. The Cayley-Menger 
determinant corresponding to N particles is given by 


X(Cl 2 , • • • jClJV) C 23 , ■ ■ ■ , C{N- 1 )n) 


= det 


0 

s 12 

Sl3 •• 

• Sin 

S 12 

0 

S23 ' ' 

■ S2N 

Sl3 

S 23 

0 •• 

■ S 3 n 

SlN 

S2N 

S3N ■ ' 

■ 0 

1 

1 

1 • • 

• 1 


(24) 


where s a p - C),:r 

In the following example we restrict ourselves to one dimension since the resulting ex¬ 
pressions are short and easy to manipulate, although this example can be readily extended 
to any dimension. It is easy to see that in one dimension the number of independent coor¬ 
dinates are N — 1 and for N > 2 the number of interatomic distances exceeds the number 
of independent coordinates. Therefore, for simplicity, consider as before a system consisting 
of three particles interacting in one dimension. The standard pair potential representation 
for this system, which is defined for all £ 12 , C 13 and C 23 is given by 


Vint(Cl2, Cl3; C 23 ) — 012(Cl2) + 013(Cl3) + 023(C23)- (25) 
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We noted earlier that the distances between particles are geometrically constrained by the 
requirement that one of the distance is equal to the sum of the other two. In spite of 
this constraint, V; nt is defined for all values of (Ci2 5 Ci3? C23)- This clearly shows that the 
pair potential is already an extension. Since the calculation gets unwieldy, let us again 
consider the special case where the particles are arranged to satisfy x\ < x 2 < x 3 , for 
which ri 3 = r\ 2 + r 2 3 - Using (l22jh the internal force, /j nt , evaluated at this configuration, is 
decomposed as 


f[ nt (r 12 ,r 13 ,r 23 ) 


dVint _ dtp 12 #13 

dx 1 dx 1 dx 1 

= #12(^12) + 0 i 3 Oi 3 ) 
=: /12 + fi3- 


(26) 


We now construct an alternate extension to the standard pair potential representation given 
in (I25|) . This is done through the Cayley-Menger determinant corresponding to a cluster of 
three points, which follows from (|24li as 

x(Ci2, C13, C23) = (C12 — C13 — C23)(C23 — C12 — C13) 

x (C13 — C23 — C12)(C12 + C13 + C23)• 


Since the Cayley-Menger determinant is related to the area formed by the three particles, 
and the three particles are restricted to be in one-dimension, it follows that 


X(ri 2 ,r 13 ,r 23 ) = 0. (27) 

Using the identity in (1271) . an alternate extension 1# is constructed: 

Thrt(Cl2, C13, C23) = Vint(Cl2, C13, C23) + X(Cl2, Cl 3 , (23)- (28) 

Note that is indeed an extension because from (12711 it is clear that is equal to Vi nt 
at every point on the shape space of the system and it is continuously differentiable because 
x(Ci2, Ci3; C23) , being a polynomial, is infinitely differentiable. Let us now see how the internal 
force, /{ nt , for the special configuration considered in this example, is decomposed using the 
new extension: 
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( 29 ) 



dV int dVint d'X 


dx i dx i dx i 



{f 12 - 8 ri 2 r 23 (ri 2 + r 23 )) + (/i 3 + 8 ri 2 r 23 (ri 2 + r 23 )) 

: fi2 + fi3, 


where in the above equation s = (V 12 , r 13 , r 23 ) is a point in the shape space S. It is clear 
from ( 1 ^ 1 ) and that the central-force decomposition is not the same for the two repre¬ 
sentations, i.e., f \2 7 ^ fi 2 and /i 3 ^ /i 3 , however the force on particle 1 , f[ nt , is the same in 
both cases as expected. 

C. Equation of Motion and the stress tensor 

The equation of motion and the stress tensor for multi-body potentials has been exten¬ 
sively studied in the authors’ previous work-22. We now present those parts of the derivation 
which are necessary to derive the energy equation in Section III D1 The equation of motion 
from continuum mechanics is given by ^ 2 



(30) 


where a is the Cauchy stress tensor and b is the body force held. Using Liouville’s equation 


and substituting in the definitions for mass density and the velocity fields defined in (J7D and 
respectively, into (l30j) . it can be shown that the stress tensor and the body force held 
must satisfy 



a 


^2 (tUV^Vint I x a = x) 


a 



(31) 


a 


where 



( 32 ) 
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is the velocity of particle a relative to the pointwise velocity field. It is natural to associate 
V ex t with the body force field b in f[3Tjl . We therefore define b(x,t ) as 

b(x, t) := - (WV Xa Vext |*a = *)- (33) 

a 

Substituting (1331) into (13TT) . we have 

cliv,,, a = - ^ m a div* (« el ® v r f)W \ x a = x) 

a 

-J2(WV Xa V int \x a = x). (34) 

a 

From (1541) . we see that the pointwise stress tensor has two contributions: 

cr(x,t) = cr k (x,t) + cr v (x,t), (35) 

where cr k and cr v are, respectively, the kinetic and potential parts of the pointwise stress. 
The kinetic part is given by 

cr k (x,t) = -^ 7n a (« el <g> v r f)W \ x a = x) . (36) 

a 

It is evident that the kinetic part of the stress tensor is symmetric. The kinetic stress reflects 
the momentum flux associated with the vibrational kinetic energy portion of the internal 
energy. 

Continuing with ([341) . the potential part of the stress must satisfy the following differential 
equation: 

diva, cr v (x, t) = y (W f" lt | x a = x) , (37) 

a 

Equation (137)) needs to be solved in order to obtain an explicit form for er v . In the original 
paper of Irving and Kirkwood^, a solution to (1371) was obtained for the special case of pair 
potential interactions by applying a Taylor expansion to the Dirac delta distribution. In 
contrast, Noll showed that a closed-form solution for <r v can be obtained by recasting the 
right-hand side in a different form and applying a lemma proved in Ref-. We proceed with 
Noll’s approach, except we place no restriction on the nature of the interatomic potential 
energy V int . 
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Derivation of the pointwise stress tensor 


We substitute the force decomposition given in (l22j) corresponding to a continuously 
differentiable extension into the potential part of the pointwise stress tensor in (157|) to 
obtain 

di v x cr v (x,t) = \x a = x). (38) 

q,/3 

<*±P 

On using the identity 

(fapW I x a = x) = / {fapW I x a = x, Xp = y) dy, (39) 

J R3 

equation (1351) takes the form 

diva, cr v (x, t) — / (Wf a/3 \ x a = x,xp = y) dy. (40) 

J R3 a 

a,p 

We now note the following lemma due to Noll, which will be used to obtain a closed-form 
solution to the output fields derived in the Irving-Kirkwood procedure. 

Lemma 1 Let f(v,w) be a tensor-valued function of two vectors v and w, which satisfies 
the following three conditions: 

1. f(v,w) is defined for all v and w and is continuously differentiable . 

2. There exists a 5 > 0, such that the auxiliary function g(v,w), defined through 

g(v, w) := f(v, it»)||u|| 3 +, 5 ||u;|| 3+<5 , (41) 

and its gradients V v g and V w g are bounded. 

3. f(v,w) is antisymmetric, i.e., 

f(v,w) = -f{w,v). (42) 


Under the above conditions, the following equation holds 



f{x,y)dy = 



f{x + SZ, X 


Us=0 


(1 — s)z) ds 


zdz. (43) 


15 




a 



FIG. 1. A schematic diagram helping to explain the vectors appearing in the pointwise potential 
stress expression in (|H]1 . The bond a-j3 is defined by the vector z. When s = 0, atom a is located 
at point x, and when s = 1, atom (3 is located at x. 

It is clear that being anti-symmetric, the integrand in the right-hand side of (HU1) satisfies all 
the necessary conditions for the application of Lemma CD Conditions (1) and (2) are satisfied 
through the regularity conditions on W. Therefore, using Lemma [H we have 


cr v (x, t) — - f I {-fapw I x a = x + sz, Xp = X - (1 - s)z) ds® z dz. 

^ a ,0 ^ R3 5=0 

a^/3 

(44) 


The expression for the potential part of the pointwise stress tensor in ()44l) is a gen¬ 
eral result applicable to all interatomic potentials. We make some important observations 
regarding this expressions below: 

1. The expression for <r v given in (144|) has an easy interpretation. <r v at a point x is 
the superposition of the expectation values of the forces in all possible bonds passing 
through x. The variable z selects a bond length and direction and the variable s slides 
the bond through x from end to end (see Fig. CO). 

2. er v is symmetric. This is clear because / aj g is parallel to z and z ® z is symmetric. 
Since the kinetic part of the stress in (FUl is also symmetric, we can conclude that the 
pointwise stress tensor is symmetric for all interatomic potentials. 

3. Since er v depends on the nature of the force decomposition and different extensions of 
a given potential energy can result in different force decompositions, we conclude that 
the pointwise stress tensor is non-unique for all interatomic potentials (including the 
pair potential). 
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The non-uniqueness of the pointwise stress tensor also plays an important role in the energy 
equation, derived in the next section, since the stress appears in it. 


D. Equation of energy balance 


The energy balance equation from continuum mechanics is given by 


dpe ,. . , 

— + div x (q - crv + pev) = 0, 


(45) 


where e is the specific internal energy and q is the heat flux vector. We saw in the previous 
section that the Irving-Kirkwood procedure extended to general interatomic potentials yields 
various possible definitions for the stress tensor. In a similar vein, we hope to use this 
extended procedure to derive possible definitions for the heat flux vector for arbitrary multi¬ 
body potentials. Before that, let us look at the definition for the heat flux vector given by 
the original Irving-Kirkwood procedure for the case of a pair potential. The heat flux vector 
in this case is decomposed as 


q := q k + q T + q v , (46) 

where 

Qk ■= \^ m a{\\ v a l f vV *W | x a = x), (47) 

a 

Qt ■■= ^ J2(vfV a W I x a = x), (48) 

a 

and 


TT-T jZ-fiapf ( ( — — - v ) w I x a = x + sz,xp = x - (1 - s)z\ ds dz , 

^ a/3 I ' Z 11 J s=0 \ V ^ J / 

q^/3 

(49) 

represent the kinetic part, transport part, and the potential part of the heat flux vector 
respectively. It was shown by Noll 3 that if the heat flux vector is defined according to (1461) . 
then along with the definition for the specific internal energy given in (lT2j) . and Lemma [Q 
the energy balance equation (145)) is identically satisfied. We can now try to extend this 
procedure to arbitrary multi-body potentials by defining a potential energy extension and 
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repeat the steps given in Ref2. But before we can do so, we must grapple with the ambiguity 
that arises in the definition for the potential part of the specific energy, e v , given in (Thfll . 
As mentioned at the end of Section III Al in order to define e v for multi-body potentials we 
must give a precise definition for the energy of each particle V a as done in dTIl) for a pair 
potential. Even for the case of identical particles, it is not a priori clear how to distribute the 
energy between the particles for a multi-body potential. It is clear from (|48|) that this results 
in an ambiguous definition for q x which depends on the definition for V a . Moreover, one 
would expect that the definitions for the pointwise fields should be invariant with respect to 
addition of any constant to the potential energy. It is clear that all the definitions discussed 
so far satisfy this invariance except for (1481) . Therefore a question that naturally arises is, 
whether the decomposition of energy is necessary to derive the energy balance equation. 
An alternate approach which we think is more reasonable conies from a paper by Murdoch” 
in his spatial averaging procedure. Here we adapt this approach to the Irving-Kirkwood 
procedure. 

An alternate derivation of the energy balance equation 

The alternate derivation for the energy balance equation in this section leads to an ex¬ 
pression for the heat flux vector which does not contain the transport part. Moreover, this 
derivation applies to any multi-body potential with a continuously differentiable extension. 
Under this alternate derivation we have the following input and output fields: 


Input Output 


' 

p 


( 


b 

V 

> ->■ < 

a 

Ck 





k Q = <lk + q v , 


We consider the terms in (1451) . beginning with pe k defined in (fT3l) . For simplicity, we assume 
Vext = 0. We have 



where E k := pe k. Using Liouvillc’s equation given in (J3]) , we obtain 
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dE k 1 


dt 




m a < tuj 


Y, ( ~ v f> ’)!*« = * 


p 


m fi 


]_ _ 1 _ 

_ 9 (ll w a|| 2 Va • Va. a W | = ®) - - ( |K \\ 2 fT ' V„ a PU \ X a = x) , 


= : <7i + <? 2 , 


(52) 


where we have used the identities (115 all and (115bl) . Now note that the term ||u Q ,|| 2 u Q , can be 
written as 


v a \\ 2 v a = |K el || 2 < el + 2« el ® 0+ 


rel ,ov „.rel\ 


2 rel 


V pud + V V 


(53) 


Consider qi, the hrst term of (1521) . Using (153)) and the definitions for q k , cr k and E k given 
in (147)) . (1361) and (fl3l) respectively, qi can be expressed as 

Qi = ~ div £C (qr k - cr k n + E k v)- 

i||n || 2 div^ (v r fW \ x a = x) 


= - div £C (qr k - er k v + E k v), 

since m a (v™ l W | x a — x) — 0. Now consider q 2 , the second term of 

by parts, and using the regularity conditions on W, q 2 takes the form 

q 2 = v a ■ ff : W | x a = x). 

a 

Using dMD and ([55]) . (1721) becomes 

dE k 

dt 

= - di v x (q k - rr k v + E k v) + • f^W \ x a = x) 

a 

= -div ;v dqi: - <r t t> + Ei,v) + ^{u™ 1 ■ f„‘W x,, = x) 

Y.(K' W ** = *> 

_ a 

We know from the momentum balance equation (see (1871) ) that 

div* cr, = i x a = x). 


(54) 
Integrating 

(55) 


V. 


(56) 


(57) 
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Using (1571) . together with the identity 


diVa;(T6) = div(T) ■ b + T . V x b, 


where T and b are continuously differentiable tensor and vector-valued functions of x re¬ 
spectively, and noting that a = cr k + er v , ([56]) can be rewritten as 


dE k 

dt 


di v x {q k - crv + E k v ) 


+ E« el • f™W \ x a = x) — cr : V x v. 

a 


(58) 


Now, consider the middle term on the right-hand side of (158|) which is given by 


*«=*>■ («*) 

a 

Substituting the force decomposition given in (122|) corresponding to a continuously differen¬ 
tiable extension, into (]59|) . we obtain 


where 


<?3 = ’ fafiW \x a = x) 

a,f3 

<*¥=& 

= 'V] [ « el • fafiW \x a = X, Xp 

a^/3 


y)dy 


= lds(x,y) + g AS (x,y)] dy, 


(60) 


9s(x,y) = 


\ ^2{{fap ■ < el + //3a ■ vf)W \ x a = x,xp = y), (61) 

ce,/3 

a^/3 


9As{x,y) = 

\ ^2((fap ■ < el - fp a vf)W \ x a = x,xp = y). (62) 

a, 0 
a^/3 

It is easy to check that g A s(x,y) = —g A s(y, x ), he., it is antisymmetric with respect to its 
arguments. Moreover, the second integrand on the right-hand side of (1601) satisfies all the 
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necessary conditions for the applications of Lemma [0 Using Lemma [H we can express the 
second integral in (|60|) as 


0as(®, y) dy = ~ cliv x q v 


(63) 


where 


*/ 

Z , o •fzeK 3 Js=0 


a„8 

af=P 


fa.fi 


V a + Vp 


v I W | x a = x + sz, xp = x — (1 — s)z ) ds dz. 


(64) 


Substituting (1631) into (l60lh and noting that 

I 9s(x, y)dy = - y^ifap ■ (v a - vp)W \x a = x), 

J R3 z 


a, P 


we obtain 


dt 


= '■ 9s(x,t) 


= - div x [q k + q v - (tv + E k v] - a : V x v + g s (x, t). 


(65) 


( 66 ) 


Now, recall the energy equation of continuum thermodynamics in (145|) . Subtracting 
form (1451) . we obtain 
d(pe. 


dt 


= - div 3 ,(q - q k - q v + pe v v) + a : X7 x v - g s (x, t). 


(67) 


The following step is a crucial part of our derivation. Note that in contrast to the original 
Irving-Kirkwood derivation, the transport part of the heat flux, q x, does not appear here. 
We can therefore identify the heat flux vector q with q k + q v , i.e., 


Thus, (1671) reduces to 


which implies that 


d(pe v ) 

dt 


q ■= qk + q v - 


= - di v x (pe v v) + cr : V x v - g s (x, t), 


( 68 ) 


dp de v . _ 

+ P~gf = U di y x {pv) - pVe v • v 

+ cr : V x v - gs(x). 


(69) 
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Using the equation of continuity, (l69]i simplifies to 

= er : X7 x v -g s (x,t), (70) 

which implies that 

pe v = cr :V x v - g s (x,t). (71) 

It is clear from (17TT) that we now have a new definition for the specific internal energy (similar 
to the one obtained by Murdoch- 7 in the Murdoch procedure) given by 

r f i 

€y(x,t)= / -{a : X7 X V -g s (x,t)) dt + c. (72) 

Jo P 

This definition does not require a decomposition of the total energy to individual atoms, 
i.e., it is independent of a particular choice for V a , contrary to what is observed in the 
original Irving-Kirkwood procedure and its generalization to multi-body potentials found 
in the literature (see Ref22— “— ). 

In summary, we obtained new definitions for e v and q , which are quite different from 
those obtained in the Irving-Kirkwood procedure. We believe that the new definitions for 
e v and q v given in (1721) and (1681) respectively are more physically reasonable as compared 
to those given in (fT4l) and (146]) due to the following features which are not observed in the 
Irving-Kirkwood procedure or any of its previous generalizations to multi-body potentials: 

1. The definitions for q and e v given in (j68j) and (1721) depend on the derivative of the 
potential thus making them invariant with respect to changes in the potential energy 
by a constant. This is a rather natural thing to expect. 

2. The heat flux vector obtained in the alternate derivation does not have transport 
part. This suggests that we look for numerical experiments which yield a non-trivial 
transport part using the original Irving-Kirkwood procedure. Most of the numerical 
experiments found in the literature, which study the energy balance equation obtained 
though the Irving-Kirkwood procedure, lump the transport part into either the kinetic 
or potential parts of of the heat flux vector and do not observe it separately. Hence, 
there has been no extensive numerical study of the role of this term. If indeed the 
expression for the transport part of the heat flux vector found in the Irving-Kirkwood 
procedure always has a negligible contribution to the heat flux vector, then its existence 


P 


aT + Vt ” 


V 
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can be questioned. Preliminary numerical simulations we conducted to explore this 
(which are not reported here) always yielded a negligible transport part. 

3. From Section Ul Gl and Section Hi Dl it follows that the only ambiguity in the expressions 
obtained through this modified derivation is the non-uniqueness of the pointwise stress 
tensor, which is directly related to the force decomposition. It was shown in Ref 2 ' that 
this non-uniqueness vanishes in the thermodynamic limit. 

III. EXPRESSION FOR MD SIMULATION 

In the previous sections we saw that various pointwise fields can be obtained through the 
Irving-Kirkwood procedure. As noted in Section II, the pointwise held are not continuum 
fields. We identify the continuum fields with macroscopic observables obtained in a two-step 
process: 

1 . A pointwise held is obtained as a statistical mechanics phase average. 

2 . A macroscopic held is obtained as a spatial average over the pointwise held. 

We have seen that the pointwise helds obtained in the hrst step are defined as phase averages 
with respect to a probability density function. Typically a molecular dynamics (MD) sim¬ 
ulation is purely deterministic in nature, meaning that at a given instant in time, we have 
a complete microscopic description of the system. Due to this knowledge, the probability 
density function introduced in the Irving-Kirkwood procedure reduces to a Dirac delta dis¬ 
tribution supported on the point in the phase space corresponding to the state of the system. 
If (a? MD (t), v MD (f)) denotes the evolution of an MD simulation, then the probability density 
function 1 D md corresponding to an MD simulation is given by 

W UT> (x,v,t) = n*<*. - “MW®. - <“(<))• (73) 

a 

Therefore, in an MD setting, the pointwise helds obtained in step 1 are localized to the 
particle positions. Next, we spatially average these helds with respect to a normalized 
weighting function that has compact support, thus obtaining expressions for the continuum 
holds that can be numerically evaluated using the data generated in a MD simulation. 
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Spatial averaging 


A macroscopic quantity is by necessity an average over some spatial region surround¬ 
ing the continuum point where it is nominally defined. Thus, if f(x,t',W ) is an Irving- 
Kirkwood- pointwise held, such as density, stress or internal energy, the corresponding 
macroscopic held f w (x,t ) is given by 

f w (x,t)=[ w(y — x)f(y, t] W) dy, (74) 

Jr 3 

where w(r) is a suitable weighting function. 

It is important to note that due to the linearity of the phase averaging in the Irving- 
Kirkwood procedure, the averaged macroscopic function f w (x,t) satishes the same balance 
equations as does the pointwise measure f(x,t). 


Weighting function 

The weighting function w(r) is a real-valued function^ with units of volume -1 which 
satishes the normalization condition 


w(r)dr = 1. 


(75) 


This condition ensures that the correct macroscopic held is obtained when the pointwise 
held is uniform. For a spherically-symmetric function, w(r) = w(r), where r = ||r||. The 
normalization condition in this case is 


w(r) = 


/ w(r)4nr 2 dr = 1. 

Jo 

The simplest choice for w(r) is a spherically-symmetric uniform function over a specihed 
radius r w , given by 

1 /K, if r < r w , 

0 otherwise, 

where V w = ~7ris the volume of the sphere. This function is discontinuous at r = r w . 
If this is a concern, a “mollifying function”- that smoothly takes w(r) to zero at r w over 
some desired range can be added (see (HOOD ). Other possible choices include for example 
Gaussian functions^, or spline function used in meshless methods^ (see Ref 27 for details). 


(76) 
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Many physical interpretations can be given to the weighting function. See Ref 7 for further 
details. 

One possible interpretation for a positive-valued w with compact support (as described 
above) can be related to the physical nature of the experimental probe measuring the con¬ 
tinuum fields. In this case, the size of the compact support represents the length scale over 
which continuum fields are being measured. An alternative approach described by Murdoch 
and Bedeaux- is based on the requirement that “repeated spatial averaging should produce 
nothing new”. In other words, spatially averaging a quantity that was already spatially av¬ 
eraged should give the same average. This leads to a definite form for the weighting function 
that also takes on negative values. 

It is straightforward to see that substituting the expression given in (1731) for the prob¬ 
ability density function into ()74j) and performing the spatial averaging defined there using 
any weighting function discussed above, we obtain expressions for continuum fields that can 
be numerically evaluated using the data generated from an MD simulation. For example, 
let us look at the mass density field given in (J7|) and repeated here with W = W MD : 


p(x,t) = J^m a (IF MD | x a = x). 

a 


(77) 


Spatially averaging this distribution with respect to the weighting function results in 


p w (x, t) = ^ m a w(x™ D - x), 

a 


(78) 


where p w denotes the continuum mass density held obtained by the spatial averaging of the 
pointwise held p with respect to the weighting function w. Similarly, all other continuum 
definitions for an MD simulation are obtained from their probabilistic versions. Following 
is a catalog of definitions for continuum helds that can be evaluated in any MD simulation: 
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p w {x,t) = m a v a w(x ] ^ D - x), 

Q 

Pw(x, t) 

V w {x,t) = - - -—, 

Pw{X,t) 

t) G"w, k(®5 T CwA X , ^)j 

e w (x, t) = e Wjk (x, t) + e WjV (x, t ), 
q w (x, t ) = q w>k (x, t ) + qr«, jV (a5, t), 

cr w , k (x, t) =- ^2 m «« el ® < e V(*a D - *), < el = - v , 


p w e w ^{x,t) = - ^m a ||'y^ ID || 2 M;(a;“ D - ®), 

a 

qw,k( x * t) = 1 5Z ™«IK"ll 2 < ll "(^r _ *). 


( 79 ) 

( 80 ) 

( 81 ) 

( 82 ) 

( 83 ) 

( 84 ) 

( 85 ) 

( 86 ) 




, t) = o / w(y-x) f (~f a pW MD I = y + sz,xp = y- (1 - s)z) ds®zdzdy, 
^ „ ft J s=0 


a,/3 


£w,v(x, t) / — (<T-u) . V Q Wt s(«C; t)) dt, g Wj s(x, t) — ^ ^ f, 

Jo P , J 


( 87 ) 




q w A x ^) = \y2 w (y 

Z «/R 3 


a (3 

a^/j 



v a +v,8 


a,P 


v ) H ,md | aj„ = y + sz, a,, = y - (1 - s)z 


( 88 ) 


( 89 ) 


The expressions for cr w v and g WiV given in (157|) and flB5j) respectively, can be further simplified 
by the following simple change of variables: 


y + sz = u, y-(l-s)z = v, 

which implies that 

z = u — V, y — {1 — s)w + sv. 
The Jacobian of the transformation is 


J = det 


V u z V„z 
j_ V„y_L 


= det 


/ -/ 

(l-s)J s/ 


= 1. 


( 90 ) 


(91) 


( 92 ) 
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Using (j90j) . (j9Tj) and (1921) . q W)W simplifies to 



OC+fi 


where 


b(x-,u,v)= / u>((l — s)u + sv — x) ds, 

J s =0 


is called the bond function. Similarly, cr WtV simplifies to 



a,(3 


Hence, equations fl75|) - fl55|l . (155)1 . (155)1 and (15511 . are the most simplified form of continuum 
fields that can be evaluated in an MD simulation. These expressions constitute a generaliza¬ 


tion to many-body potentials of expressions originally given by Hardy- for the special case 
of pair potentials (noting that our expressions for the internal energy density and heat flux 
vector which have modified forms). Other commonly-used definitions like the virial stress^ 
and the Tsai traction^ can be obtained as limiting cases of the above relations (see Ref 2 ' 


for details). We therefore see that our formulation is unified in the sense that it shows how 
all of these relations are related and provides a single framework that encompasses them all. 

IV. NUMERICAL EXPERIMENTS 

It was mentioned at the end of Section III Dl that the definitions obtained in the new 
derivation are quite different and some qualitative differences were identified. We now try 
to see how the definitions vary quantitatively. Various stress tensors obtained through this 
unified framework were studied in Ref 27 by the authors. In this section, we describe molecular 
dynamics simulations that are carried out to compare the following two quantities: 



a 


(95) 



o P 


1 


((Tw : V x v w - g WiS (x,t)) dt. 


(96) 
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Here ej^ v is the local spatial average of the potential part of the specific internal energy in 
the original Irving-Kirkwood formulation. e w>v is the corresponding expression in the new 
formulation taken from (I88H . The constant in (1951) is chosen in order to compare the above 
two equations from a fixed datum. To explore the role of the two parts of the integrand in 
(ESI) ) we define 


4,vO^) 

4,v0m) 


so that e WiV = e } w v + e 2 wy . 



• V x^w) dt, 



-g w ,s(x,t) dt , 
p 


(97) 

(98) 


Interatomic potential 


Since the unified framework applies to arbitrary multi-body potentials, it would have been 
ideal to choose a multi-body potential for our numerical experiments. Unfortunately, since 
as mentioned in the introduction, there is no rigorous way to distribute the total energy 
among particles, the expression for e^ v in (1951) is ill-defined. Thus the only possibility 
for comparison is in the special case of pair potential interactions in a system of identical 
particles. In this case due to symmetry, it is reasonable to divide the energy equally among 
the particles (see footnote^), thus arriving at a plausible definition for e^ v . We will see that 
even in this case the expressions in (|95l) and (|96l) are not equivalent. In more general cases, 
we argue that only the new expression in (1961) is well-defined. 

Our system consists of particles arranged in a face-centered cubic crystal, stacked together 
in the form of 15 x 15 x 15 unit cells, interacting through a modified Lennard-Jones type 
potential. The Lennard-Jones parameter, e and cr are arbitrarily set to 1. The potential has 
the following form: 


</>(r) = 4 


y>12 /v>C 

The above equation has been rendered dimensionless by expressing lengths in units of cr and 
energy in units of e. As seen in the above equation, the standard Lennard-Jones potential 
is modified by the addition of a quadratic term. This is done to ensure that </>(r cut ) = 
0 / ( r cut) = 0, where r cut = 2.5, denotes the cutoff radius for the potential. We refer to this as 
the “Modified Lennard-Jones potential”. The ground state of this potential is an fee crystal 


- 0.0078r 2 + 0.0651. 


(99) 
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FIG. 2. Plot showing the average specific potential energy and the average specific kinetic energy. 
Since this is a constant energy simulation, their sum (black solid line ) is always constant. 


with a lattice constant of a = 1.556. Tims, the dimension of the sample at ground state 
is given by its length l = 15 x 1.556 = 23.34. We use the velocity Verlet time integration 
algorithm to evolve the system. 

The weighting function, w(r), is chosen to be a constant with a suitable mollifying 
function,- 


w(r) 


c if r < R — 5 

\c [l — cos ( ; yy r 7r)] if R — 5 < r < R , 


( 100 ) 


0 


otherwise 


where 5 = 0.12, R is the radius of the sphere which forms the compact support and c = 
1 — (5 / R) 3 + 3 (5/R) 2 — 1.5(5/R). The constant c is chosen to normalize the weighting 


function. 


Experiment 1 

We begin with a constant energy molecular dynamics simulation with periodic boundary 
conditions. The atoms in the system are randomly perturbed to induce a temperature of 
T = 0.145 in Lennard-Jones units after equilibration. Fig. [2] shows the evolution of the 
average specific potential energy, i.e., the total potential energy divided by the mass (= N 
in our case), and average specific kinetic energy, which add up to a constant specific internal 
energy. (The average specific internal energy is constant since this is a constant energy 
simulation.) 

Now, suppose we are interested in evaluating the continuum quantities given in (195]) and 
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(a) 


(b) 


FIG. 3. Evolution of specific internal energy for a constant energy simulation with periodic bound¬ 
ary conditions, (a) Plot showing the evolution of the potential part, e WiV , and the kinetic part, 
e «j,k) °f the specific internal energy. The total specific internal energy (shown in black solid line) is 
not strictly constant, (b) Plot comparing the evolution of the potential part of the specific internal 
energy with its analogue, e)j) v , in the original Irving-Kirkwood procedure. 


at the center of our system. To do this, we choose R = 0.4/, which corresponds to a 
length of 6 unit cells, for the radius of the compact support of the weighting function. 


Fig. 3(a) shows a plot of the potential and kinetic part of the specific energy at the center 
of the sample calculated using the weighting function given in (HOOIh It is clear from the plot 
that the total specific energy at the center of the sample has some oscillations up to about 
200 time steps before these oscillations become negligible. This reminds us that although we 
are performing a constant energy simulation, the specific internal energy at a point need not 


be constant. Fig. 3(b) compares the different expressions for the potential part of the specific 
internal energy given in (1951) and (1961) . The plots for e WtV and v are in good agreement. It 
is also clear from the plot that the contribution of e l w v to e WtY is negligible. 


Experiment 2 

In this experiment, we continue with the same system with periodic boundary conditions. 
We begin with the particles at their equilibrium positions, with temperature T = 0. We allow 
the system to evolve for the first 1000 time steps during which we observe small fluctuations 
due to numerical noise. We then instantaneously increase the temperature to T = 0.2 using 
a simple velocity rescaling thermostat and maintain this temperature for the rest of the 
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(a) 


(b) 


FIG. 4. Evolution of specific internal energy for a constant temperature (applied after first 1000 
time steps) simulation with periodic boundary conditions, (a) Plot showing the evolution of the 
potential part, e W}V , and the kinetic part, e w ^, of the specific internal energy (b) Plot comparing 
the evolution of the potential part of the specific internal energy with its analogue, e^ v , in the 
original Irving-Kirkwood procedure. 


simulation. Again, we are interested in studying the continuum fields at the center of the 
sample. We choose the radius of the compact support R = 0.4/ for the weighting function, 


as in the previous experiment. Fig. 4(a) shows the plot of e W}V , e w ^, and the total specific 
internal energy at the center of the sample for this case. Fig. 4(b) compares e WjV with e* v . 
It is clear from this plot that both e W)V and e^ v are in good agreement with each other. It 
is interesting to see that the contribution to e W:V due to el vy remains negligible even after 
increasing the temperature of the system. 


Experiment 3 


The setup is similar to the previous experiment except that now we do not apply periodic 
boundary conditions. This means that the sample is free to expand once the temperature 


of the sample is increased. The plots shown in Fig. [5] correspond to R — 0.4/. Fig. 5(a 


shows the plot of e WtV , e W) k and the total specific energy at the center of the sample, and 
Fig. 5(b) compares e WiV with v . In this case we see from Fig. 5(b) that unlike the previous 
two experiments there is a small negative contribution to e w , v due to v — although the 
magnitude is still very small. Note that the longer wavelength oscillations in energy in these 
plots correspond to the pulsing of the sample as it expands and contracts about its mean 
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(a) (b) 

FIG. 5. Evolution of specific internal energy for a constant temperature (applied after first 1000 
time steps) simulation without periodic boundary conditions, using an averaging domain of radius 
R = OAl. (a) Plot showing the evolution of the potential part, e WjV , and the kinetic part, e w ^, of 
the specific internal energy, (b) Plot comparing the evolution of the potential part of the specific 
internal energy with its analogue, e^ v , in the original Irving-Kirkwood procedure. 


thermally-expanded size. 

It is also interesting to see how these continuum fields change as the averaging domain 
size is decreased. As mentioned previously, as the weighting function tends to the Dirac 
delta distribution, we expect the continuum fields to also become localized with increasing 
magnitude. To see if this is indeed the case, we decrease the size of the averaging domain 
by choosing R = 0.1/, which corresponds to a length of 1.5 unit cells. Fig. O shows the 


plots for this case. Comparing Fig. 6(a) with Fig. 5(a), we see that e w ^ remains the same, 
while e w>v for the smaller averaging domain is about four times larger than its value for the 


larger domain. Similarly comparing Fig. 6(b) with Fig. 5(b), we see that for the smaller 
domain, e w , v is greater than four times e^ v , whereas they were equal for the larger domain. 
Moreover, the contribution due to v to the specific internal energy is larger than that due 
to e^ v . This clearly shows that for small averaging domains, the two definitions given by 


(1951) and fl96l) are quite different in nature. Based on our observations, we can conclude that 
e WjV tends to localize more strongly with the averaging domain size than does e^ v . 

To rationalize these results, let us consider what happens to these definitions as the 
weighting function tends to a delta distribution. In this case, e^ v (see (J95])) localizes to 
the particle positions. The same also happens to all terms in the definition of e WtV (see 
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(a) (b) 


FIG. 6. Evolution of specific internal energy for a constant temperature (applied after first 1000 
time steps) simulation without periodic boundary conditions, using an averaging domain of radius 
R = 0.1/. (a) Plot showing the evolution of the potential part, e WjV , and the kinetic part, e w ^, of 
the specific internal energy, (b) Plot comparing the evolution of the potential part of the specific 
internal energy with its analogue, ej]j v , in the original Irving-Kirkwood procedure. 

(|96ll ) except for cr w . This term localizes to the lines joining the particles rather than onto 
the particles themselves. This observation provides a qualitative explanation, although not 
complete, for the different behavior of the two definitions shown above. 

V. SUMMARY 

In this paper, we present a two-step unified framework for the evaluation of continuum 
field expressions from molecular simulations: (1) pointwise continuum fields are obtained 
using a generalization of the Irving-Kirkwood procedure to arbitrary multi-body potentials, 
and (2) spatial averaging is applied to obtain the corresponding macroscopic quantities. The 
latter lead to expressions suitable for computation in molecular dynamics simulations. It is 
shown that the important commonly-used microscopic definitions for continuum fields are 
recovered in this process as special cases. 

In generalizing the Irving-Kirkwood to arbitrary many-body potentials we have had to 
address two ambiguities inherent in the original procedure which lead to non-uniqueness 
in the resulting expressions. The first ambiguity arises due to the non-uniqueness of the 
partitioning of the force on an atom as a sum of central forces, which is directly related 
to the non-uniqueness of the potential energy representation in terms of distances between 
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particles. This is in turn related to the shape space of the system. The conclusion is that the 
pointwise stress tensor is not unique, however we show in Ref 2 - that the macroscopic stress 
obtained via spatial averaging becomes unique as the spatial averaging domain is taken to 
infinity. 

The second ambiguity in the original Irving-Kirkwood procedure arises due to the ar¬ 
bitrary decomposition of energy between particles. We show that this decomposition can 
be completely avoided through an alternative derivation for the energy balance equation. 
This leads to new definitions for the specific internal energy and the heat flux vector. In 
particular, the resulting potential part of the specific internal energy does not depend on 
the arbitrary partitioning of the potential energy to individual particles, and the resulting 
heat flux vector does not contain the “transport part” which is not invariant with respect 
to the addition of a constant to the potential energy function. 

The new definition for the specific internal energy is compared with the original Irving- 
Kirkwood definition through a series of numerical experiments. Although our expression 
applies to arbitrary many-body potentials, we have chosen to perform the comparisons for 
the special case of a system of identical particles interacting through a pair potential since 
this is the only case where the original Irving-Kirkwood internal energy density is well- 
defined. This is due to the ambiguity in the decomposition of energy between the particles 
in the original Irving-Kirkwood derivation (and existing extensions to the approach). In 
our numerical experiments, we observe that both definitions agree for weighting functions 
with support given by a ball of radius 0.6 unit cells and larger. However, as the weighting 
function tends to a delta distribution, the two definitions scale differently. A qualitative 
theoretical explanation for this difference is given based on the limiting behavior of the two 
definitions as the averaging domain tends to a point. 

It would also be of interest to compare the new expression for the heat flux vector to the 
original Irving-Kirkwood expression. In order to do this, one has to study the transport 
part of the heat flux. To our knowledge, this has not been done in the past because the 
transport part is normally lumped into either the kinetic or potential parts of the heat flux 
and not observed separately. Our new derivation shows that its existence is directly related 
to the arbitrariness in the energy decomposition. Preliminary numerical studies of some 
systems (not reported here) suggest that the transport part of the heat flux vector, which is 
absent in the new derivation, tends to be negligible. Further work is necessary to determine 
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whether this is a general result. 


ACKNOWLEDGMENTS 

This research was partly supported through the National Science Foundation (NSF) under 
Grant No. DMS-0757355 and the Air Force Office of Scientific Research (AFOSR) under 
Grant No. FA9550-09-1-0157. The views and conclusions contained herein are those of 
the authors and should not be interpreted as necessarily representing the official policies or 
endorsements, either expressed or implied, of the AFOSR or the U.S. Government. 


REFERENCES 

1 J. H. Irving and G. Kirkwood, The Journal of Chemical Physics 18, 817 (1950). 

2 The derivation is non-rigorous in the sense that expressing the stress tensor as a series 
expansion is only possible when the probability density function, which is used in the 
derivation, is an analytic function of the spatial variables (see Ref~). 

3 W. Noll, Journal of Rational Mechanics and Analysis 4, 627 (1955). 

4 R. Lehoucq and A. Von Lilienfeld-Toal, Journal of Elasticity 100, 5 (2010). 

5 R. J. Hardy, The Journal of Chemical Physics 76, 622 (1982). 

6 A. I. Murdoch and D. Bedeaux, International Journal of Engineering Science 31, 1345 
(1993). 

7 A. I. Murdoch and D. Bedeaux, Proceedings: Mathematical and Physical Sciences 445, 
157 (1994). 

8 A. I. Murdoch, Journal of Elasticity 71, 105 (2003). 

9 A. I. Murdoch, Journal of Elasticity 88, 113 (2007). 

10 Although the Hardy and Murdoch procedures seem similar, they are notably different in 
the derivation of the energy balance equation, and the resulting expressions for energy 
density and the heat flux vector. 

n A. Cauchy, “Exercises de mathematique,” (Chez de Bure Freres, Paris, 1928) Chap. Sur 
l’eequilibre et le mouvement d’un systeme du points materiels sollicites par des forces 
d’attraction ou de repulsion mutuclle, pp. 227-252. 


35 


12 A. Cauchy, “Exercises de mathematique,” (Chez de Bure Freres, Paris, 1928) Chap. De 
la pression ou tension dans un systeme de points materiels, pp. 253-277. 

13 R. Clausius, Philosophical Magazine 40, 122 (1870). 

14 See Ref^. for a more detailed historical review. 

15 D. H. Tsai, The Journal of Chemical Physics 70, 1375 (1979). 

16 E. Wajnryb, A. R. Altenberger, and J. S. Daliler, The Journal of Chemical Physics 103, 
9782 (1995). 

17 J. Cormier, J. M. Rickman, and T. J. Delph, Journal of Applied Physics 89, 99 (2001). 

18 J. A. Zimmerman, E. B. Webb, J. J. Hoyt, R. E. Jones, P. A. Klein, and D. J. Bammann, 
Modelling and Simulation in Materials Science and Engineering 12 , S319 (2004). 

19 H. Heinz, W. Paul, and K. Binder, Physical Review E 72, 066704 (2005). 

20 T. J. Delph, Modelling and Simulation in Materials Science and Engineering 13, 585 
(2005). 

21 S. Morante, G. C. Rossi, and M. Testa, The Journal of Chemical Physics 125, 034101 
(2006). 

22 Y. P. Chen, The Journal of Chemical Physics 124, 054113 (2006). 

23 K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, The Journal of Chemical Physics 
130, 204106 (2009). 

24 K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, Physical Revew E 80, 047702 
(2009). 

25 A. P. Thompson, S. J. Plimpton, and W. Mattson, The Journal of Chemical Physics 131, 
154107 (2009). 

26 G. C. Rossi and M. Testa, The Journal of Chemical Physics 132, 074902 (2010). 

27 N. C. Admal and E. B. Tadmor, Journal of Elasticity 100, 63 (2010). 

28 J. A. Zimmerman, E. B. Webb, and S. C. Seel, Mathematics and Mechanics of Solids 13, 
221 (2008). 

29 J. Zhang and B. D. Todd, Physical Review E 69, 031111 (2004). 

30 G. Marechal and J.-P. Ryckaert, Chemical Physics Letters 101, 548 (1983). 

31 D. Torii, T. Nakano, and T. Ohara, The Journal of Chemical Physics 128, 044504 (2008). 

32 The symmetry argument for equally dividing the energy among particles only holds for 
identical atoms interacting via a pair potential. It is lost for multi-species systems and for 
all higher-order potentials. For example, for a three-body potential (even in the case of a 


36 


single-species system), the symmetry between atoms is lost in all clusters of three particles 
which do not form an equilateral triangle. The same reasoning applies for potentials of 
higher order. 

33 The usual convention is to represent the phase space via positions and momenta of the 
particles. For convenience, in this section, we instead use positions and velocities. 

34 The assumption that the probability density function exists and it is continuously differen¬ 
tiable can be considerably weakened by viewing IT as a generalized function/distribution 
in the sense of Schwartz. For the sake of brevity we do not take this approach, however, we 
later use a generalized function/distribution as a candidate for W to arrive at expressions 
for continuum fields that can be used in a molecular dynamics simulation. See Section IIIII 
for details. 

35 If G is a second-order tensor or higher, then the dot product indicates tensor operating 
on a vector. Note that in (fT5]h in the interest of brevity, we are breaking our notation of 
denoting a second-order tensor operating on a vector by juxtaposition. 

36 Note that this assumption may fail in systems undergoing first-order magnetic or electronic 
phase transformations. 

37 E. B. Tadmor and R. E. Miller, Modeling Materials: Continuum, Atomistic and Multiscale 
Techniques (Cambridge University Press, Cambridge, 2011). 

38 M. Daw and M. Baskes, Phys. Rev. B 29, 6443 (1984). 

39 H. Stillinger and A. Weber, Physical Review B 31, 5262 (1985). 

40 J. Tersoff, Phys. Rev. B 38, 9902 (1988). 

41 We thank Ryan Elliott for suggesting this line of thinking. 

12 L. E. Malvern, Introduction to the mechanics of a continuous medium (Prentice-Hall, 
Upper Saddle River, 1969). 

43 The expression in Noll’s paper appears transposed relative to (l43]h This is because the 
gradient and divergence operations used by Noll are the transpose of our definitions. 

44 It was mentioned in Ref 27 that the weighting function is a positive-valued function based on 
our interpretation of it being related to the nature of the experimental probe. We thank 
1. Murdoch for directing us to an alternate interpretation which allows the weighting 
function to take on negative values. See Ref 7 for details. 

45 T. Belytchko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, Computer Methods in 
Applied Mechanics and Engineering 139, 3 (1996). 


37 


